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ABSTRACT 

i.  ■  4  .  •  •  - 

We  discuss  a  method  for  suppressing  the  oscillations  of  a  linear  system 
subject  to  an  external  periodic  disturbance  of  fixed,  but  unknown,  period. 

The  method  entails  augmentation  of  the  original  plant  with  a  compensator  and 
parameter  identifier.  The  near  equilibrium  dynamics  of  the  resulting  system 
are  analyzed  and  shown  to  be  related  to  a  linear  delay  equation  with  infinite 
delay  and  periodic  coefficients.  A  corresponding  Floquet  theory  is 
indicated.  A  FORTRAN  program  approximately  realizing  the  period  identifier  is 
included  and  numerical  results  obtained  with  this  program  are  graphically 
displayed  and  analyzed. 


MOS:  93B30,  B40,  C40,  45A05,  D05. 

Key  Words:  Adaptive  Control,  Decoupling,  Delay  Equations 
Work  Unit  Number  5  -  Optimization  and  Large  Scale  Systems 
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SIGNIFICANCE  AND  EXPLANATION 


For  a  wide  variety  of  systems,  including  sighting  devices,  weapons, 
machine  tool  arms,  etc.,  operation  under  conditions  which  involve  significant 
oscillatory  disturbances  is  necessary.  Often  it  is  desirable  to  dynamically 
decouple  the  system  from  the  disturbances  by  means  of  the  intervention  of 
active  control.  In  many  cases  this  must  be  done  without  a  prior  knowledge  of 
the  period  (equivalently,  the  frequency)  of  the  incoming  disturbance.  In  this 
paper  we  propose  a  method  for  such  vibration  suppression  using  a  compensator 
and  frequency/period  identifier.  The  stability  of  the  resulting  complex  is 
analyzed  and  numerical  studies  are  presented  to  indicate  the  potential 
effectiveness  of  the  procedure. 
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FFFCUENCY/PERIOD  ESTIMATION  AND  ADAPTIVE 
REJECTION  OF  PERIODIC  DISTURBANCES 

D.  l.  Russell* 


0.  INTRODUCTION 


In  a  wide  variety  of  applications  one  encounters  a  system  of  the  form 

x-Ax+Cu  +  v,  (x-  ll")  (0.1) 

wherein  x  Is  the  n-dimenslonal  state  vector,  u  Is  the  m-dimensional  control  vector 
and  v  Is  a  periodic  n-dlmenslonal  vector  disturbance  function  with  least  positive 
period  T: 

v(t)  -  v(t  +  T)  .  (0.2) 

In  many  cases  x  »  Ax  by  Itself  represents  the  dynamics  of  an  elastic  system,  the 
disturbance  v  arises  from  the  environment  in  which  the  elastic  system  is  placed,  and  the 
control  u  is  used  to  mitigate  the  effects  of  this  disturbance.  Examples  include 
sighting  devices  (cameras,  telescopes,  etc.),  weapons,  and  machine  tool  arms,  operated 
under  conditions  which  involve  significant  oscillatory  disturbances,  such  as  would  be  the 
case  for  a  telescope  operated  from  an  aircraft,  e.g..  Another  important  application 
arises  in  connection  with  the  measurement  and  active  suppression  of  aerodynamic  flutter  in 
aircraft  wings,  tail  structures,  etc. 

The  approach  taken  in  this  paper  is  to  suppose  that  v(t)  can  be  adequately  modelled 
by 


v(t)  -  B*(t) ,  B  n  *  2r , 
where  z(t)  satisfies  a  linear  system 

Z  “  Ft  , 


(0.3) 


(0.4) 
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the  kj  being  positive  integers.  These  need  not  necessarily  be  1,2,  ...,r;  in  some  cases 
it  is  known,  e.g.,  that  only  odd  order  harmonics  occur  so  that  we  would  use  k.j  =*  1, 
k2  =  3,  ...,  kr  =  2r  -1. 

Assuming  F  known,  and  this  will  bring  us  to  the  subject  of  frequency  estimation 
later  on,  we  can  construct  a  compensator 

y  =  Sx  +  Fy  (0.7) 

where  y  is  the  2r-dimensional  compensator  state,  and  consider  the  combined  system 
x  =  Ax  +  Bz  +  Cu 

y  =  Sx  +  Fy  (0.8) 

z  =  Fz  . 

We  will  suppose  that  the  range  of  C  includes  the  range  of  B.  This  means  that,  in 
principle,  one  could  solve 

Cu  =  -Bz  (0.9) 

and  cancel  the  effect  of  the  disturbance  altogether.  For  a  telescope  operated  from  a 

moving  vehicle,  neglecting  translational  motion  and  considering  only  the  angular 
displacements,  this  would  be  the  case  if  the  controls,  acting  through  the  mounting,  have 
both  azimuth  and  elevation  correctional  capability.  In  practice  the  direct  cancellation 
(0.9)  is  rarely  feasible  due  to  noise,  measurement  delays,  limited  measurement  capability. 


Let  4  be  any  nonsingular  2c  *  2c  matrix  which  commutes  with  Fi  In  most  cases  we 
would  use  the  identity  matrix.  We  may  then  find  an  m  *  2r  matrix  L  such  that 

CL  -  -B*"1  .  (0.10) 

Assuming  additionally  that  (A,C)  is  stabilizable,  let  K  be  an  m  x  n  matrix  such 
that  A  +  CJC  is  a  stability  matrix  and  let  u  be  generated  by  the  feedback  control  law 

u  -  Kx  +  Ly  .  (0.11) 


Using  this  in  (0.S)  we  have 


we  will  see  in  Section  1  that  it  is  possible  to  select  S  in  such  a  way  that  the 
control  law  (0.11)  dynamically  decouples  the  plant  state  x  from  the  periodic 
disturbance  v(t)  •  Bz(t). 

The  foregoing  scheme,  to  be  developed  more  fully  in  the  next  section,  clearly  amounts 
to  the  construction  of  a  reduced  order  observer  for  the  disturbance  state  s(t)  (see 
[10])  and  assumes  that  the  plant  state  x(t)  is  completely  accessible.  If  this  is  not 
the  case,  dynamic  decoupling  is  probably  best  realized  with  the  construction  of  a  full 
n  +  2r  dimensional  state  observer.  Assuming  an  observation 

w  -  HqX  +  H1z  (0.13) 


available  such  that  the  pair 


(H0,H,)  [  A  +  CK  B  j 

\  o  r) 


(0.14) 


-3- 


the  estimator  system 


\  =  (A  +  CK)5  +  Lg  ( Hg x  -  HgC)  +  Lg  ( H  ^  z  -  ) 

r.  =  Fz  +  L^HgX  -  HgC)  +  )  . 

Then,  choosing  u  such  that 


Cu  =  -Bf 

and  letting  e  =  x  -  r,  F  =  z  -  F  we  find  that 


and  we  conclude,  since  (0.15)  is  a  stbility  matrix,  that 

lim  e(t)  «  lim  f(t)  =  0  . 
t-«»  t'*40 


Since,  with  (0.19),  (0.16),  (0.17)  become 
x  =  (A  +  CK)x  +  Bf 
z  *  Fz 

we  conclude  that 


lim  x(t)  »  0 

t» 


(0.18) 


(0.19) 


(0.20) 

(0.21) 


and  thus  x(t)  is  decoupled  from  z.  This  is  a  standard  procedure,  such  as  described  in 
[10],  for  example. 

Whether  decompling  is  carried  out  as  in  Section  1  or  as  above,  it  is  clear  that  the 
estimator  system  requires  knowledge  of  the  matrix  (0.5)  and  hence  the  parameter  a  =  2*/T 
in  (0.6).  When  the  period,  T,  and  hence  a,  is  unknown  it  is  necessary  to  adjoin  a 
parameter  estimator  to  supply  the  system  with  an  estimate  for  T.  Such  a  parameter 
estimator  is  described  in  Section  2.  Stability  considerations  in  connection  with  the 
period  estimator  lead  to  examination  of  a  related  functional  equation  of  retarded  type  in 
Section  3.  A  numerical  realization  of  the  estimator  of  Section  2  is  developed  in  Section 
4  and  examples  of  its  use  are  presented  in  Section  5. 


1.  COMPENSATOR  DESIGN  FOR  A  KNOWN  DISTURBANCE  FREQUENCY 

If  the  period,  T,  or,  equivalently,  the  frequency  v  »  1/T,  of  the  disturbance  v 
is  known,  then  we  may  assume  that  F  is  known  and  the  only  problem  in  constructing  the 
compensator  (0.7)  is  the  selection  of  the  2r  x  n  matrix  S.  Let  us  note  that  the  matrix 


equation 


fA  t  CK  -B»  jfOj  _  f01p.  +  fB)  ,  o 
S  f  ft  *  o 


is  clearly  valid  whatever  S  may  be.  This  means  that  if  we  define  n  by 


we  shall  have 


n -  rA  + «  -B*  .  (i. 

n  S  F  n 

as  is  easily  checked.  If  the  matrix  in  (1.3)  is  a  stability  matrix,  then  x(t)  ”  £(t) 


will  have  the  property 


lia  «x(t)»  -  0 
t-">» 


so  that  the  periodic  disturbance  v(t)  “  8z(t)  has  only  a  transient  effect  on  x(t)/  the 
range  of  the  transfer  function  matrix  from  z  to  includes  only  vectors  of  the  form 

f®].  Thus  the  plant  state  vector  x  is  dynamically  decoupled  from  z.  If  the  matrix  in 
(1.3)  is  not  a  stability  matrix  no  such  inferences  are  valid,  our  proof  that  S  can  be 
selected  so  as  to  satisfy  this  stability  requirement  begins  with 


THEOREM  1  Let  F  be  antihermltlan,  (as  in  (0.5)),  so  that 

F*  -  -F  . 

Then  the  n  x  m  linear  matrix  equation 


(A  +  CK)PQ  -  PqF  -  B®  =  0 


has  a  unique  n  x  2r  solution  P0.  If  the  pair  (P0,F)  is  observable,  then  the  2r  *  n 


matrix  S  can  be  choaen  in  such  a  way  that 


M 


fA  +  CK 


(1.5) 


is  a  stability  matrix. 


Proof  Since  F*  =  -F  implies  that  F  has  only  purely  imaginary  eigenvalues,  the 
existence  of  a  unique  solution  PQ  of  (1.4)  is  assured  by  a  familiar  theorem  in  matrix 
theory  (see,  e.y.,  (2]).  An  easy  application  of  the  implicit  function  theorem  then  shows 
that  the  cubic  matrix  equation 

(A  +  CK)P  -  PF  -  B*"1  +  EPP*P  =  0(P,£)  =■  0  (1.6) 

has  a  unique  solution  P  *  P(e)  defined  for  small  e  ♦  0  with 

lim  P(e)  =  P  (1.7) 

£♦0 

Setting 


S  =  S(e)  =  -ep(e)* 
we  note  that  M  in  (1.5)  is  similar  to 


1  -P(e)  A  +  CK  -B*-1  I  P(e ) 

B‘*>  *  f  S  «._  ^-ep(e)*  F  "  "  -  1 


2r 


0  I 


2r 


(-A  +  CK  +  ep(e )p(e  )*  P(P(E),E)  1 

-ep(e)»  F  -  ep(e)*p(e) 


/■A  +  CK  +  ep(s)P(e)*  0  \ 

1  -ep (€)*  F  -  ep(e)P(e)» 


(1.8) 


Since  K  has  been  chosen  so  that  A  +  CK  is  a  stability  matrix, 

M  (E)  5  A  +  CK  +  EP(E)P(£)» 
n 

is  an  n  x  n  stability  matrix  for  sufficiently  small  E  >  0.  From  the  antihermitian 
property  of  F  we  can  see  that 

O  -  £P(E)*P(£)1*  I.  +  I,  ( F  -  £P( E ) *P( E ) 1  +  2EP(£)*P(E)  =  0. 


Applying  a  well  known  modification  of  Liapounov's  Theorem  (See,  e.g.,  (81)  we  conclude 


that 

M2r ( e )  ~  F  -  cP(e )*P(e ) 

is  a  stability  matrix  for  £  >  0  if  (P(e),F)  is  observable.  Since  we  have  assumed 
(Pa,F)  observable  and  (1.7)  is  true,  (P(e),F)  is  observable  for  e  >  0  sufficiently 
small  and  M2r(£)  thus  a  stability  matrix  for  these  values  of  E,  at  least.  Since 

M(e)  ia  lower  block  triangular  with  blocks  Mn(e),  M^fE),  *tB  atabmty#  ®nd  hence 
that  of  M  «  M(£)  in  (1.5)  is  assured  with  the  choice  (1.8)  for  S  «  S(e)  for 
sufficiently  small  e  >  0. 

It  will  be  noted  that  the  choice  of  the  feedback  matrix  K  is  import  .  n  at  least 
two  ways.  Improvement  of  the  convergence  of  *  V  to  z,  i.e.,  reduction  of  t>  - 
transient  effect  of  the  disturbance  v  =  Bz,  dictates  choosing  e  larger  to  improve  the 
stability  properties  of  F  -  ep(e)*p(e).  But,  since  A  +  CK  +  £P(e)P(e)*  suffers, 
stability-wise,  as  £  is  increased,  K  must  be  UBed  to  offset  this  effect.  In  Section  3 
we  will  find  even  further  considerations  to  take  into  account  in  the  selection  of  e 
and  K. 

Since  P(e)  satisfies  a  cubic  equation,  which  may  entail  some  difficulty  of 
solution,  the  following  corollary  is  useful  in  applications. 


(1.9) 


(1.10) 


Proof  With  the  indicated  choice  of  S  the  matrix  M( E )  is  similar  to 

(x  -p0\/a  +  ck  -b^Wi  p\ 

\  o  I  /  \  -EP0*  F  /  \  0  I  / 


-7- 


A  +  CK  +  epQpg*  (A  +  CK)PQ-  PoP  -  B®  +  EP0Po*Po' 

-ep  *  F  -  ep  *P„ 

0  0  0  i 


=  cf.  (1.4)  = 


A  +  CK  +  EP0P0*  ^p0p0*p0 

-£P0*  F  -  £P0*P J  • 


=  ,V2 


for  e  >  0.  with  PQ(g)  =  uPq  •  the  matrix  (1.11)  becomes 


A  +  CK  +  P0(U)P0(U)*  -  P0(u)P0(U)*P0(u) 


-gP0 ( u ) * 


F  -  P0(g)*P0(u) 


which  is  similar  to 


A  +  CK  +  PQ (u )P0 (U  )*  P0(u)P0(u)*P0(u) 


-P0(U)* 


A ..  (U )  P3(U) 

-pQ(g)*  F1  ( u ) 


F  -  pQ(g)*P0(u) 


The  corresponding  lower  triangular  matrix 


-P0(g)*  F^u) 


is  a  stability  matrix,  using  essentially  the  same  argument  as  in  Theorem  1,  provided 


g  >  0  is  sufficiently  small.  Consider  the  equation 


1  ( g )*  -PQ ( u  Aj  Q  f  \  [P  Rl[  A^g) 


0  F 1  (  g )  ^  R*  T  J  ^  R*  T /\-PQ  (U  ) *  F ^  ( g  ) 


f 1  °  \  -  ° 

I  o  2P0(u)*P0(g)  I 


Solving  this,  we  find  that  T  =  I2r 


A ^  ( g ) *Q  -  p 0 ( g ) R*  +  QA^g)  -  RPQ ( g ) *  +  1  =  0 


A. (u)*R  -  P„(g)  +  RF , ( g ) 


0 


< 


I 

p 

I 


For  small  u  (equivalently,  small  e  )  the  eigenvalues  of  A^Cu)  and  -F^fu)  are 
uniformly  separated  and  solution  of  (1.15)  shows  that 

R  -  C>f  IP0(u)ll  -  0<U> 

and  then  a  similar  analysis  of  the  first  equation  shows  that 


0  “  0Q(u)  +  £)(u  ) 


where 


VlO'OoOO  +  CgdOA^u)  +  In  -  0  . 

Thus  Qg(U),  and  hence  p,  remains  bounded  for  u  >  0  small.  Since  (1.13)  is 
satisfied,  using  the  matrix  of  (1.12)  instead,  we  have 


(a, (u )*  >p0(u)Ue  r\  /e  rV a,(u)  p3(m 

P3<U)*  t)  +  1r«  -P  (U)*  ^  (U 


1  7rcVo*' 


o  /i  up0*  /  \  v0*°  J2  VP 


+  I  QP„P„*P„P„ 
1  2  *00  0  0 


*0  0  0  | 

0/ 


From  this  it  is  easy  to  see  that  the  matrix  on  the  right  hand  side  is  nonpositive  for 
small  v  >  0  and  the  result  follows  by  the  familiar  Liapounov  theorem,  provided  that 
whenever 


fih 


VU) 

-P0(U) 


* 


P3(U) 

F I  (U  ) 


-9- 


the  quadratic  form 


cannot  vanish  on  any  interval  of  positive  length.  For  small  u  >  0  this  question  reduces 
very  quickly  to  the  observability  of  the  pair  (Pg,F),  which  has  already  been  assumed. 
This  completes  the  proof. 


2.  PERIOD  ESTIMATION  AND  A  RELATED  STABILITY  PROBLEM 

Whether  the  submatrix  S  in  the  definition  (1.5)  of  M  is  selected  as  in  Theorem  1 
or  as  in  Corollary  2,  or  by  some  other  procedure,  it  is  clear  that  the  overall  matrix  M 
will  depend  on  the  period,  T,  of  the  disturbance  v  so  that,  supposing  now  that  the 
design  parameter  e  has  been  fixed,  we  have 


M  »  M(a) 


th  +  CK  -B*(a)-1  ■) 
S(a)  F(n) 


(2.1) 


where  <*  *  — (if  (0.6)). 

It  would  be  possible  to  estimate  a  directly  using  various  well  known  parameter 
estimation  procedures  ([6],  [9],).  However,  in  these  procedures  one  tends  to  encounter 
either  instability  or  slow  convergence,  or  other  difficulties.  For  example,  the  model 
reference  algorithm  of  [6]  cannot  be  applied  because,  in  the  complete  system  (0.12)  the 
portion  z  *  Fz  of  that  system  is  not  controllable  with  respect  to  u. 

We  have  elected  to  use  a  very  simple  procedure  to  estimate  the  period,  T, 
directly.  Assuming  that  an  output,  or  observation 

di(t)  -  H,x(t)  +  H2y(t)  -  'W^tP  S  Hw(t) ,  (2.2) 

where  (cf.  (0.3),  (0.12),  (1.5)). 

w  »  M(«)w  +  R  («  -  (p))  (2.3) 

is  available,  from  the  assumed  stability  property  of  the  matrix  M(n)  it  follows  that  a 
periodic  disturbance  input  v  will  result  in  an  output  w(t)  which,  except  for  transient 
behavior,  is  also  periodic  with  the  same  period.  It  therefore  makes  sense,  in  the 
continuous  framework  which  we  use  here  for  analysis,  to  consider  the  cost  functional, for 
v  >  0, 

CQ( T,t)  *  f  eY*S  fc'fw(s)  -  w(s  -  T)Ww(s)  -  w(s  -  T)1ds 
0 

*  f  eY^S  t'fw(s)  -  w(s  -  Tg 1 *H*h( w( s)  -  w(s  -  T)lds 
0 


-11- 


and  select,  as  our  estimate  for  the  period  T  at  time  t,  that  value  T(t)  which 
minimizes  C(T,t)  within  a  given  range  T^  <  T  <  Tj  •  (The  range  must  be  restricted  in 
order  to  avoid  the  trivial  period  T  =  0  and  multiples  of  the  minimum  period  of  the 
disturbance.)  Then 


a(t) 


2w 

T(t) 


(2.5) 


is  the  estimate  at  time  t  for  a  in  (0.6),  (2.3).  A  numerical  procedure  approximating 
this  optimization  process  is  described  in  Section  4  and  is  used  to  obtain  the 
computational  results  of  Section  5.  There  it  will  be  seen  that  certain  steps  do  have  to 
be  taken  in  order  to  ensure  the  stability  of  the  combined  control/estimation  system.  Our 
purpose  here  is  to  provide  a  framework  for  the  stability  analysis  by  developing  a 
linearized  variational  equation  for  that  system  about  the  nominal  time  trajectory  in  the 
case  where  the  true  period,  which  we  will  call  Tg,  lies  in  the  interior  of  the  interval 
T1  <  T  <  T2. 

For  our  analysis  of  the  combined  use  of  (2.3)  and  (2.4)  we  will  consider,  instead 
of  Cg (T, t ) ,  as  given  by  (2.4),  the  cost 


C(T,t)  «  f  eY(8't)fw(s)  -  w(s  -  T) w( s>  -  w(s  -  T) )ds  , 
0 


(2.6) 


wherein  we  assume  that  the  trajectory  w(s)  is  defined  in  the  indefinite  past, 
justification  for  this  lies  in  the  fact  that  if  (2.3)  and  (2.6)  together  yield  a 
process,  the  difference  between  the  use  of  (2.6)  and  (2.4)  will  be  transient. 

To  carry  out  this  program  we  begin  by  supposing  that  when  the  correct  value 
is  used  in  (2.3)  the  steady  state  Tg  -  periodic  solution  resulting  from  the  T0 
periodic  input  v(t)  is  wg(t)  and  u>g(t)  “  Hwg(t).  Since  our  estimate,  a(t), 
from  we  suppose  that  the  actual  solution  of  (2.3)  which  we  obtain  is  w(t) 

Pit)  -  °(t  -  Tq) 

W(t)  =  Wg ( t)  +  ^W(t)  , 


The 

stable 


will  vary 
Thus 

(2.7) 

(2.8) 
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«{t)  -  «0  +  An(t) 


(2.9) 


T(t)  -  TQ  +  4T(t)  .  (2.10) 

A  necessary  condition  in  order  that  T(t)  should  minimize  C(T,t)  is  obtained  by 
differentiating  C(T,t)  with  respect  to  T  and  setting  the  result  at  T  -  T(t)  equal  to 
zero.  Thus  (cf.  (2.2),  (2.3),  (2.5)) 

0  -  1  .(?!*)■  I 

2  »T  -  T{ t) 

■  (  eY*s  t'fw(s)  -  w( 8  -  T(t) )1*H*Hw(s  -  T(t ) )ds 

—m 

rn  I  eY*®  t^w(8)  -  w(s  -  T(t)  ))*H*f A(«(s  -  T(t)  ))wf  s  -  T(t)1 

+  pfs  -  Tttj’Uds  .  (2.11) 

Noting  (2.5),  we  see  that  (2.11)  is,  implicitly,  an  equation  for  T(t)  which  is  coupled 
with  the  system  (2.3)  satisfied  by  w(t).  The  resulting  coupled  system  is  clearly  a 
nonlinear  functional  equation  of  delay  type.  He  are  concerned  with  the  (at  least  local 
with  respect  to  o0  and  wQ ( t ) )  existence,  uniqueness  and  asymptotic  stability  of 
solutions . 

LEMMA  3.  for  fixed  t  and  a  trajectory  w(s)  <  s  <  t,  for  (2.3),  corresponding  to 
a  continuous  T0  -  periodic  P(s),  -  *  <a  <  t,  the  equation  (2.11)  is  solvable  for 
T(t)  near  Tq  ^f_ 

-  fw(t)  -  w(t  -  T_ ) )*H*w( t  -  T  ) 

0  0 

V  (  a-f  )  •  I 

+  '  e  w(s)*H*w(s  -  TQ)ds 

+  y  f  eY^8  t^fw(s)  -  w(s  -  Tq  )  1  *H*w(  a  -  Tg)ds 


+  a(  eY  <8-t>(w<s)  -  v(»  -  Tq ) )*H*Hw( s  -  TQ)dslAT(t) 

+  I  eY'*  tVw(»)  -  w(s  -  Tq )  1*H*Hw( ■  -  Tq)(Jb  +  ....  (2.17) 

This  expression  no  longer  depends  on  w(s  -  TQ),  hence  we  may  relax  the  requirement 
Pec1  since  flee0  can  be  uniformly  approximated  in  the  C°  norm  by  fl  e  C1. 
Examination  of  the  remainder  term  shows  that  the  same  argument  applies  there  and  we 
conclude  that  (2.16)  is  indeed  valid  to  first  order  in  AT(t).  The  firat  statement  of  our 
lenmia  then  follows  imnediately  from  the  implicit  function  theorem. 

The  last  statement  follows  from  the  property 

w0(s)  -  wQ(s  -  Tq)  =0,  <  s  <  t  (2.18) 

together  with  (2.3),  (2.7),  (2.17).,  which  enable  one  to  make  the  first  term  in  (2.16) 
arbitrarily  close  to  the  left  hand  side  of  (2.15)  and  the  second  term  arbitrarily  close  to 
zero.  This  completes  the  proof  of  Lemma  3. 

The  linearization  with  respect  to  AT,  Aw  is  obtained  by  using  (2.18)  in  (2.16). 
Because  wQ( s)  -  wQ( a  -  TQ)  3  0  and  because  only  zero  order  terms  are  retained  aa 
coefficients  of  the  firat  order  term  AT(t),  the  result  is 

rr  eY("_t>  w0<s)*H*Hw0(s  -  Tq )ds 1 At( t ) 

+  f  eY*8  tVAw(s)  -  Aw(s  -  Tq))*H*Wq(s  -  Tg)ds  “  0  , 

•  I 

and  using  (2.18),  wQ(s  -  TQ )  =  w(s),  and  the  assumption  (2.15)  we  have 

-f  eY*B  ^(Awts)  -  Aw(s  -  Tq )  1*H*HWq( s)ds 

AT(t)  -  — — - - -  ,  (2.19) 

(  eY(8_t)  Wg(s) *H*HWq( s)ds 

which  is  a  delay  type  functional  equation  relating  AT(t)  and  the  time  history  of  Aw. 
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Then  from  (2.3)  we  have,  to  first  order. 

Aw  +  Wg  ■  M(ag)Wg  +  6  ♦  M(aQ)Aw  +  (aQ)Aa)w0  +  ... 

and,  since  wQ  »  M*a0*w0  +  g  and  (2.5)  applies,  we  have  as  our  linearized  system  (2.19) 
and 

Aw(t)  -  M( aQ )Aw(t )  -  [ZL-  |£  (a0)wQ(t)jAT(t)  .  (2.20) 

T0 

A  single  equation  may  be  obtained  by  noting  that  for  t  >  TQ 


Aw(s)  -  Aw(s  -  Tq) 


M(aQ) (Aw(s)  -  Aw( s  -  Tq ) )  -  ||  (a0)w0(s)AT(s) 

T0 


It  <VV"  -  VAT<*  -  To) 

To 


M(a  )(Aw(s)  -  Aw(s  — T  ) )  -  (a  )w  (s)(AT(s)  -  AT(s  -  T  )  , 


*0''  m  2  3a  "*0'w0' 
T0 


where  we  have  used  w0 ( s )  «  wQ ( s  -  Tq ) •  Then 


Aw(t)  -  Aw(t  -  Tq) 


t  M(a  )(t-o), 

~  /  e  (a0)w0(o)(AT(a)  -  AT(a  -  TQ))ds 


and  thus 


AT(  t) 


/  Y ( S~t )  •  • 

j  eT  Wp(s)*H*Wp(s)ds 


rC  T(s-t)  ,*  *  3M  .  *  ««*0  >*<«-«> 

j  e  J  wo(0)  to  ‘“o’  6 


(AT(O)  -  AT( 0  -  TQ ) ] do  H*HwQ(s)ds  . 
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Letting 


W0<Y,t>  -  — 


Y(s-t)  • 


w0(e)*H*Hw0(e)de 


(2.21) 


we  have 


t  *  aw  .  t  .  .  M(a  )  (e-o) 

AT(t)  -  W  (Y.t)  /  w  (a)  ,4  (a  )  J  eY1  'e  H*Hw  (s)ds 


V°'  JS  Ka0 


*  (AT(O)  -  AT(0  -  Tq  )  ) do  , 


which  has  the  form 


AT(t)  -  /  W1(Y,t,O,M(O0))(AT(0  -  Tq ) ] do 


(2.22) 


W^Y.^O.Mlc^))  -  W0<Y,t>W0(of  («0  )*  /  eY(8_t)e  ° 


M(cr.)  (a-o) 


H*Hwq (a )ds  (2.23) 


LEMMA  4.  w1(T»t,o,M(a0))  is  periodic  in  t  and  o  with  period  T0,  in  the  sense 

W^Y.t.o^OQ))  -  W^Y.t  +  T0,O  +  T0,M(aQ))  . 

Proof.  Using  the  formula  (2.23)  directly  we  see  that 

W^Y.t  +  T0,o  +  T0,M(a0))  -  W0(Y,t  +  T0)w0(O  +  TQ)*  §£  (a Q> * 


t>Tft  Y(s-t-T„)  M(o  )  (8-O-T  ) 

/  0  e  0  e  0  0  H«Hw0(s)ds  . 

o+T0 


From  the  T0  -  periodicity  of  w  (t)  the  same  periodicity  of  W (Y.t)  follows  easily. 


Then  with  r  «  s  -  Tn 


W,(Y,t  +  Tq.O  ♦  T0,M(a0))  -  W0(Y,t)w0(o)*  U  (Oq)1 


*  Y(r-t)  M<a0>  <t-°) 

J  e  © 


H*Hwft(r  +  T  )dr 


♦ 


and,  using  the  Tq  -  periodicity  of  we  have 

W^Y.t  +  TQ,c  +  T0,M(a0)l  -  W^Y.t.o.MInQ))  , 

as  claimed. 

The  equation  (2.22)  can  be  rewritten  as 
t 

AT(t)  «  f  wfY,t,3,M(a0)UT(C)d(j  (2.24) 

with 

wfY/t,3,M(a())’l  »  W1fY,t,c,M(a0)''/  t  -  Tq  <  a  <  t 

wfY,t,a,M(a0))  «  W^Y^.CjMtag  )1  -  W1  f  y  #t  ,o  +  T0,M(a0>),  -•  <  a  <  t  -  TQ  . 

Since,  for  a  -  t  -  TQ,  W^Y.tjd  -*■  tq,  M(oq))  »  W1  (y  ,t,t,M(aQ ) )  «  0  we  see  that 
wfY,t,3,M(aQ) )  is  continuous  as  a  function  of  a  and,  clearly, 
wf Y ,t  +  T0,c  +  T0,M(a0))  -  wfY,t,o,M(a0))  . 

We  have  shown  Y,  M(aQ)  directly  as  arguments  of  W  because  y,  e  and  K,  the  last  two 
involved  in  the  construction  of  M( aQ )  (see  (2.1))  are  the  parameters  which  we  have  to 
work  with  in  order  to  influence  W,  and  hence  the  solutions  of  (2.24).  It  is  clear  from 
(2.22)  that  W  depends  on  Tq  as  well. 
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ANALYSIS  OF  FLOQUET  TYPE  SOLUTIONS 


l 

i 


The  fact  that  the  equation  (2.24)  involves  an  infinite  time  delay  places  it  in  a 
class  of  functional  differential  equations  with  periodic  coefficients  whose  properties 
have  not  been  fully  explored.  From  the  behavior  of  solutions  of  such  equations  with 
finite  delays  ([3],  [4])  we  expect  that/  with  some  restrictions  on  the  kernel 
w[Y,t,0,M(ao)),  the  dominant  solutions  should  be  solutions  of  "Floquet  type",  i.e., 
solutions  of  the  form 

AT(t)  -  eXtP(t)  ,  (3.1) 

where  P(t)  is  a  continuous  Tq  -  periodic  functions 

P(t  +  T0)  «  P(t)  . 

The  main  point  of  this  section  is  to  indicate  that  this  is,  indeed,  the  case  for  kernals 
satisfying  a  uniform  decay  condition 

|w(y,t,c,M(a0)] j  <  c  e"c(t-0)  ,  a  <  t  , 

for  positive  C,  c. 

Before  entering  upon  the  proof  of  this,  let  us  note  some  rather  transparent  results 
which,  however  superficial,  give  some  indication  of  the  factors  which  are  likely  to  play  a 
role  in  our  analysis.  Suppose  an  inequality  (3.2)  is  satisfied  for  positive  c,  C. 
Supposing  a  solution  of  the  form  (3.1)  to  exist,  we  normalize  P(t)  so  that 

sup  | P( a > |  -  1  . 
se[t,  t  -  Tq] 

Then  we  let  t  be  such  that  |P(t)|  «  1.  Multiplying  by  a  constant,  if  need  be,  we  may 
assume  P(t)  “  1.  Then 

eXt  -  j  w(r,t,o,M(a0))eXoP(a)do  -  1 

—CD 

or 

/  w(r,t,0,M(ao))eX<o-t)P(o)do  -  1 
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V 


But 


MY,t,o,M<c.n))eX(0-t)P<°)|  <  C  e^c  +  Re(X>)(0-t) 


jVc  +  R«(X))(o-t)do  „  ,  f 


yielding  an  upper  bound  on  Re ( X ) : 


c  +  Re(X) 


>  1  “>  Re(X)  <  C  -  c  . 


Under  what  circumstances  could  a  bound  of  the  type  (3.2)  be  expected?  Recalling  that 
W1(Y/t,o,M(o0))  -  wo(Y,t)w0(a>*  H  (a0)V<o-t) 


t  y (s-o )  M(V*(8-0)  • 

«  j  e  e  H*HWg(s)ds 

a 


we  note  that  with  r  “  s  -  a, 

y(s-o)  (M<V*  +  •  t“a  (M(a  )•  +  YI)r  • 

1  e' '  e  H*Hw(r  +  o)ds  -  j  e'  1  o'  ’  H*Hw  (r  + 


Since  wQ  is  periodic,  if  the  eigenvalues  u  of  M(aQ)  satisfy 

Rely)  <  -  6  , 


for  some  6  >  0,  we  will  have,  for  some  Mg  >  0 
(m(oq)*  +  Yl)r  . 


H*HwQ(r  +  a)  I  <  MQe 


(Y-5  )r 


t-o  (M(a.)*  +  YI)r 
Ij  e 
0 


H*HwQ(r  +  o)drl 


<  M  I  e  dr  -  — j-  [e  -  1]  • 

0  1 


We  expect  WQ(t,Y)  to  be  0(y)  from  (2.21);  write 

|w0(t'-r)l  4  M1Y 

and  then,  since  Wg(o)  is  periodic, 

»w0^a,*  "5a  (a0**e^a-t>l 

<  M2Y  eY(a_t>  .  (3.4) 

Combining  this  with  (3.3),  (2.23) 

|w,(Y,t,O,M<a0))|  <  MqM2  JCy  e610-0 

giving  C  -  M„M2  ,  c  -  5  . 

A  comparable  estimate  will  then  apply  to  w(Y»t,0,M(ajj ) ) . 

From  this  we  see  that  if  we  are  to  control  the  identifier  stability  properties,  this 

must  be  done  through  Y  and  through  the  system  matrix  M(Yq),  by  choice  of  y ,  t  and 

K  (or  through  choice  of  Y>  K  ,LQ ,  L1  if  we  use  the  full  system  estimator  as  decribed  in 

Section  0.  Further,  we  see  from  (2.21)  that  W0<Y,t),  and  hence  W^y  ,t,o,M(aQ )) , 

w(Y»t,a,M(a  ))  increase  rapidly  as  the  frequency  parameter  a  -  increases,  i.e., 

T0 

as  Tq  decreases.  Thus,  to  be  able  to  reject  higher  frequencies  while  maintaining 
stability  we  must  expect  to  find  it  necessary  to  increase  the  damping  in  the  system  (2.3) 
by  use  of  higher  gains  e  and  K  (of  (1.5),  61.8)).  We  will  also  see  in  Section  5  that 
this  expectation  is  realized. 

We  proceed  now  to  state  a  theorem  to  the  effect  that  if  a  bound  of  the  form  (3.2) 
applies,  then  all  solutions  of  (2.24)  which  do  not  satisfy 

IAT(t) I  <  Be"6t,  0  <  t  <  «  ,  (3.5) 

where  B  is  positive  and  6  >  0  is  less  than  c  by  an  arbitrarily  small  amount,  must  be 


linear  combinations  of  Floquet  type  solutions 


THEOREM  5.  Consider  the  vector  functional  equation 


z(t)  =*  /  W(t,s)z(s)ds,  z  e  Rm, 


where  W(t,s)  is  a  (piecewise  continuous,  at  least)  m  x  m  matrix  function  satisfyin 


I w( t , s ) I  <  Ce-c  t_s  ,  -»  <  s  <  t  ,  ( 

W(t  +  T,  s  +  t)  =  W(t,s)  ( 

for  positive  numbers  C,  c,  T.  Then,  given  any  6  <  c,  and  any  solution  z(t)  with 
locally  square  integrable  initial  function  satisfying 
0  2 

/  e2CSlz(s ) I  ds  <  »  ,  (3 

—09  R1" 

we  can  write 


Z(t)  “  Zp(t)  +  Zg(t),  t  >  0 
where,  for  some  positive  B, 


(3.10) 


lz0(t)l  <  Be  p  .  t  >  0 

p 


(3.11) 


and  z„(t)  is  a  linear  combination  of  Floquet  type  solutions,  i.e.,  solutions  of  the  form 


C(t)  *  eXtP(t),  P(t  +  T)  »  P(t),  pec  (fO,»J,RmJ  , 


(3.12) 


or,  in  some  cases  (multiple  “eigenvalues") 


Ut)  =  eXttpP(t)  , 


(3.13) 


where  p  is  a  positive  integer  and  P(t)  is  as  in  (3.12). 

A  complete  proof  of  Theorem  5  is  beyond  the  scope  of  the  present  work  but  a  sketch  of 
the  proof  will  be  given  in  Section  6. 

From  this  result  we  see  that  whenever  an  inequality  of  the  type  (3.2)  is  valid  with 
c  >  0,  then  all  solutions  of  (2.24)  decay  at  a  uniform  exponential  rate  unless  there  are 
actually  solutions  (3.1)  of  Floquet  type  for  which  Red)  >  0.  The  question  arises,  of 


course,  as  to  how  such  Floquet  exponents  might  actually  be  computed.  It  seems  almost 


certain  that  the  most  efficient  procedure  involves  actual  solution  of  (2.24)  or  (2.19) 


(2.20),  assuming  an  adequate  approximation  procedure  is  available.  The  procedure  is 
essentially  the  same  one  as  is  used  to  compute  the  dominant  (pairs  of)  root (a)  of  an 
ordinary  polynomial. 

Returning  to  AT(t)  as  the  name  for  the  solution,  we  select  a  more  or  less  arbitrary 
initial  state  AT(t)  on  some  interval  £  — x  # 0 ]  (in  terms  of  Section  6  this  should  be  a 
z  such  that  the  residue  of  (l  -  fi(X))  1  q(X,z)  at  X  •  Xz  is  not  zero,  which  is 
generically  true).  The  resulting  solution  AT(t),  t  >  0,  is  confuted  and  we  examine 
successive  segments  of  length  Tq 

AT^( s)  -  AT(kTQ  +  s),  0  <  s  <  Tq,  k  -  0, 1,2,3, 

If  the  largest  multiplier 

X  T„ 
v  0 

Mv  -  e 

is  a  unique  real  number,  then  generically  with  respect  to  the  choice  of  initial  function 
AT(t),  t  e  [— T , 0 ] ,  we  shall  have  (using  the  least  squares  approach) 


u 


v 


lim 

k+» 


T 

j  ATk(s)ATlc_1(8)ds 


J  (AT.  (s))2ds 
0 


In  the  case  of  a  dominant  complex  conjugate  pair  the  procedure  is  only  slightly  more 
complicated.  He  solve 

ATk(s)  +  aATk_1(s)  +  BATk_2(s)  -  0 

for  a  and  B  in  the  least  squares  (least  L2  norm)  sense,  which  amounts  to 
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The  pair  Mv»  Wv  is  then  approximated  at  the  k-th  stage  by  the  roots  ,  P^  of 

u2  +  a,  p  +  8,  =  0  . 
k  k 

It  seems  likely  that  while  (2.24)  is  nicer  from  the  viewpoint  of  mathematical  simplicity, 
it  is  better  to  solve  (2.19),  (2.20)  rather  them  (2.24)  because  the  formula  for  the  kernel 
w(Y,t,<J,M(a0 )  )  in  (2.24)  is  rather  complicated. 

If  a  simulation  routine  combining  the  period  estimator,  compensator  and  a 
methematical  model  of  the  plant  to  be  controlled  is  already  in  hand,  as  was  the  case  for 
the  writer,  approximate  solutions  of  the  variational  equation  can  be  obtained  by  running 
the  simulator  with  slightly  different  initial  conditions  and  forming  the  appropriate 
difference  quotient  of  the  resulting  solutions.  This  'oes  not  test  the  validity  of  our 
derivation  of  the  variational  equation  but,  as  we  will  see  in  Section  5,  it  does  provide 
results  consistent  with  the  proposed  functional  equation  model  for  error  propagation. 


NUMERICAL  REALIZATION  OF  THE  PERIOD  ESTIMATOR 


It  x(t)  is  a  solution  of 

x(t)  “(At  at)  x(t)  +  v(t),  t  >  0  , 
and  tha  disturbance  v(t)  is  periodic  with  period  T: 

v(t)  -  v(t  +  T)  , 


then  an  observation  on  x(t) , 


<u(t)  •  H  x{ t ) 

will  tend  eiqponentially  to  a  period  observation,  i.e., 

li»  («(t)  -  w(t  ♦  T))  “  0  . 

t*** 

In  this  section  we  develop  a  numerical  procedure  for  estimation  of  T  which  is  a 
realization  of  the  continuous  procedure  described  in  Section  2.  We  will  take  u  to  be 
scalar  here  but  the  extension  to  vector  observations  is  quite  immediate. 

We  will  suppose  that  u>(t)  is  not  available  continuously.  Rather,  we  have  discrete 
samples 


“k  -  “(tk>#  *k+1  -  fck  *  h  >  °'  *  “  °'1'2 . 

For  computational  purposes  we  define  the  interpolated  observation  on  t^  <  t  <  tfc+1  by 

2(tk  +  oh)  =  n^to)  -  0“k+1  ♦  (1  "  0)«k,  0  <  o  <  1  .  (4.1) 

We  note  that  «»(tk)  “  w(tk+1 )  “  •  We  hk  “  k  «  0,1,2 .  Our  method 

for  estimating  T  is  to  form,  at  each  instant  tk,  and  for  a  range  LQ  <  l  <  L1 ,  the 
functions 


pk.t(0) 


\  "  \-t(0> 


and  determine  values  l,  ,  o,  of 

k  k 


1,  o  which  minimize 


(4.2) 


*  : 


(0)  -  l  vV  ,  ,(0)  , 

.*  J.0  k  3,* 


which  should  be  compared  with  (2.6).  The  functions  (4.2),  of 
values  nk  ”  id^,  nk+1  “  for  this  description  and  the  p 


course,  require  only  the 
k  l  admit  a  comparable 


we  see  that 


I 


s 

i 

ii 

> 

. 

i 


Each  pair  4,  o  corresponds  to  a  delay  t(4,0) 
associated  with  a  function  C^lt)  defined  for  t^  - 
4  corresponding,  when  a  -  0,  to  the  points  t  - 


we  pass  from 

t 

-  tk  -  4h 

to  t  -  t^  -  (4  - 

3Ck,i 

.  !fk 

3o 

o  - 

0 

3t 

t 

-  (tk-4h)+ 

3Ck,4  | 

3Ck 

L 

3o  ' 

o  - 

1 

3t 

-  (t  -  ( 4-1 )h)-  . 

candidate  for  the  minimizing  value  just  in  case 


■  (4  -  o)k.  Thus  < o  >  can  be 

Ljh  <  t  <  t^  -  L^h,  the  values  of 
t^  -  4h.  As  a  increases  from  0 
1)h.  Thus  we  have 


_  >  0  1  3Ck,i+1 

2  3d  I  '2  3<J  .  , 

'a  «  0  'o-1 


1  3Ck,4 


<  0  , 


i  •  G  •  f 


Sk-4  "  Pk,k-4+1  +  Pk,k-4  +  Pk-4+1,k-4  >  0  ' 


Sk-4  _  Pk,k-4  +  Pk,k-4+1  “  Pk-4,k-4-1  *  0  * 


On  the  other  hand,  the  interval  [(4  -  1)h,  4h]  is  a  candidate  for  containing  the 
minimizing  value  of  Tk  just  in  case 


1  3C*.* 

2  3a 


i  •  ©  •  / 


1  3Ck  4 
>0  —  *» 

•  2  3  a 

a  -  1  o  »  0 


<  0  , 


Sk-4+1  "  Pk,k-4+1  +  Pk,k-4  *  Pk-4+1,k-4  >  °  * 


(4.6) 


Sk-4  "  P  k,k-4+1  +  Pk,k-4  +  Pk-4+1,k-4  <  0  * 


(4.7) 


If  (4.6),  (4.7)  are  true  for  a  given  4,  we  compute  the  corresponding  o  by  setting 
(4.5)  equal  to  zero,  i.e., 

lSk-4  ~  Pk, k-4+1  *  Pk, k-4  *  Pk-4+1,k-fJ 


[Sk-4+1  “  Sk-4  '2Pk-i+1,k-i] 
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Once  the  finitely  many  possible  candidates  for  have  been  selected  by  this  process, 

is  chosen  froa  these  as  the  one  yielding  the  smallest  value  of  ^(o). 

It  is  possible  to  economize  on  memory  space  by  using  slightly  modified  quantities. 

With 


Pk,4 


■  n.  -  n, 


k-t 


Sk,i  “  ^  Y  lPk-j,4)  ' 


Pk,l,i-1  Y  Pk-j , 4pk- j ,4-1 


^  Y3(pk-j,t-1)2  ' 


updated  via 


’k+1,4  "  (nk+1  “  V  +  pk,4-1 '  1  “  1,2»  •••'  L2' 


Sk>1,4+  (Pk*1.t)  +YSk,t  '  *  ■  '•* . L2 


Pk+1,4,4-1  “  Pk+1,4pk+1,4-1  +  Y  Pk,4,4-1,4  “  2*  L2  ' 


it  may  be  seen  that  we  have 


C.  ,(o)  -  o2[s.  -  2P  ,  ,  ,  ♦  S.  .  ,1  -  2<j[s.  ,  -  P.  ,  „  J  +  S.  „  ,j 

K  f  X>  K 0  X.  Kf  ft  f  K|  l*  1  Kf  ft  Kf  ft  |  ft*l  Kf  ftvl 


.  3c,  .  .  .  *  .  * 

—  — ~J—  »  ofs  -  2P  +  S  1  -  fs  -  P  1 

2  3o  1  k,4  k, 4,4-1  k,4-V  L  k,4  k,4,4-1J 


and  this  vanishes  when 


The  other  aspects  of  the  analysis  remain  as  above.  This  procedure  is  the  one  actually  used 
in  Fortran  SUBROUTINE  PERIOD  (L2,  LI,  GAMMA,  H,  PER,  Y),  whose  listing  follows  and  which 
forms  the  basis  for  the  numerical  experiments  carried  out  in  Section  5.  Here  PER  is  the 
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returned  period  estimate  while  Y  is  the  supplied  observation  at  each  instant  tk 
other  arguments  are  parameters  whose  identification  is  clear. 


The 


SUBROUTINE  PERI0D(L2. LI. CAMMA, H,  PER, Y) 

DIMENSION  S ( 40 ) , P  <  40 ) , RHO ( 40 ) 

INTEGER  L. LI, L2. L1P1,  LM1 

L2M1  =  L2  -  1 

ETAOLD  =  ETANEW 

ETANEW  =  Y 

DO  IS  L  -  1, L2M1 

KL  =  L2  -  L  +  1 

KLM1  =  KL  -  1 

15  RHO(KL)  =  ETANEW  -  ETAOLD  «-  RHO(KLMl) 

RHO ( 1 )  =  ETANEW  -  ETAOLD 

DO  16  L  =  1,  L2 

16  S (L )  =  GAMMA*S(L)  +  (RH0(L>)**2 
DO  17  L  =  1, L2M1 

LP1  =  L+l 

17  P  <  L  )  =  GAMMA*P  <  L )  +  RHO < L > *RH0 ( LP 1 ) 

PER  =  FLOAT (LI ) *H 

SMINI  =  S<  L 1 > 

L IP  1  =  Ll  +  1 
DO  26  L  =  L IP  1 ,  L2 
)  LM1  =  L-l 

)  IF(S(LM1 ).  LT. SMINI )G0  TO  21 

)  GO  TO  22 

)  21  PER  =  FLOAT ( LM1 ) *H 

)  SMINI  =  S(LM1) 

'  22  IF(P(LM1 >  GT  S(LM1 ) >G0  TO  24 

IF ( P (LM1 ) .  GT.  S ( L ) >G0  TO  24 
QUO  =  S  <  L ) +S ( LM 1 ) —2.  *P(LM1) 

IF(QU0  LT.  000001 )G0  TO  26 
SIG  =  (S(L)-P(LMl ) )/QUO 

VALSIG  =  SIG*SIG*S(LM1 )  +  ( 1.  -SIG>*< 1  -SIG >*5(L>+2.  *SIG* 
1(1  -SIG ) *P ( LM1 ) 

IF(VALSIG  LT  SMINI )G0  TO  23 
GO  TO  24 

23  PER  =  (FLOAT (L) -SIG )*H 
SMINI  =  VALSIG 

24  IF(S(L).  LT  SMINI )G0  TO  25 
GO  TO  26 

25  PER  =  FLOAT ( L  >  *H 
SMINI  =  S ( L ) 

26  CONTINUE 
RETURN 
END 


If  there  is  some  danger  of  confusing  the  minimal  period  T  with  one  of  its 
multiples  2T,  3T,  etc.,  this  can  usually  be  overcome  by  specifying  LI,  L2  correctly. 
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i.  SOME  NUMERICAL  EXPERIENCE  WITH  PERIOD 


We  have  carried  out  extensive  computer  based  simulations  using  SUBROUTINE  PERIOD 
described  in  the  preceding  section.  Here  we  will  describe  results  obtained  in  connection 
with  the  plant 


♦  (5.x-  * 

X  X 

with 

v(t)  -  z1 (t) 


(5.1) 


(5.2) 


1%)  “  l-o  O)  ^2^  *  (5‘3) 

t  z 

The  rather  nondescript  parameters  appearing  in  (5.1)  result  from  the  fact  that  this  damped 

inertial  system  is  a  model  for  a  certain  physical  plant  of  interest.  A  compensator  was 

* 

constructed  in  the  form  (0.7)  using  S  “  -ePQ  with  P0  as  in  (1.10).  To  provide  a  sore 
or  less  standard  basis  of  comparison  the  feedback  coefficients  in  all  cases  were  chosen  to 
achieve  critical  damping  (i.e. ,  multiple  real  eigenvalues)  at  various  rates  in  the  closed 
loop  matrix 


A  +  CK 


(  °  1  ) 

V28k1  -.2+28k2J 


In  this  very  simple  example 


(5.4) 


B  -  C  - 


4  and  -L  (cf.  (0.1))  were  both  chosen  to  be  2  *  2  identity  matrices.  For  the  values 
of  ki  and  k2  which  were  used  Pq  turns  out  to  be  a  very  small  matrix  and  we  used 
e  «  1000. 

The  output  used  for  the  period  estimation  was  w(t)  “  x^(t)  +  y1 (t )  and  the  output 
shown  on  the  diagrams  is  x^(t).  It  will  be  seen  that  the  initial  estimates  for  the 


period,  T(t),  are  wildly  inaccurate  but,  in  the  cases  when  the  complete 
plant/compensator /identifier  system  is  asymptotically  stable,  the  estimate  T(t) 


converges  to  the  correct  value  Tq  (within  the  accuracy  permitted  by  the  approximations 

inherent  in  PERIOD,  as  described  in  the  preceding  section).  In  all  cases  we  selected 

* 1 ( 0 )  »  1,  z2  ( 0 )  “0,  so 

1  2 

z  (t)  -  cos  (at),  z  (t)  -  -a  sin  (at)  . 

In  Figures  1-16  the  odd  numbered  Figures  show  the  output  x'ft)  while  the  next, 
even-numbered  figure  in  each  case  shows  the  period  estimate  T(t)  for  the  same  run. 

Figures  1  through  4  correspond  to  choices  of  k^  and  k2  such  that  the  matrix  (5.4) 
has  a  double  eigenvalue  X  «  -5.  In  Figs.  1,  2  a^  =  20*  (10  Hertz)  corresponding  to 

Tq  ■  .1,  indicated  by  the  dotted  line.  Figures  3,  4  illustrate  the  corresponding 

experience  for  aQ  «  30*  (15  Hertz),  or  Tq  "  .0667.  Here  the  period  identifier  diverges 
from  the  correct  value  and,  as  seen  in  Fig.  3,  no  significant  reduction  of  the  oscillation 
of  x'(t)  is  realized.  We  believe  that  this  is  accounted  for  by  the  fact  that  the  term 
2*/(Tq)2  in  (2.21)  changes  from  200*  in  the  10  Hz  case  to  450*  in  the  15  Hz 
case . 

Figures  5-8  show  the  10  Hz  case  with  k^  and  k2  chosen  so  that  (5.4)  has  a 

double  eigenvalue  X  «  -8  (Figs.  5  and  6)  and  with  k1  and  k2  chosen  so  that  X  “  -9. 

These  cases  seem  quite  satisfactory  with  rapid  attenuation  of  the  oscillation  in  x1 (t) , 
better  in  the  second  case  than  in  the  first,  and  rapid  convergence  of  the  period  estimate 
T(t)  to  Tq  »  .1.  The  corresponding  experience  in  the  15  Hz  case  is  not  nearly  so 
satisfactory.  Figs.  9  and  10  show  the  performance  for  X  *  -8.5  while  Figs.  11  and  12 
show  X  =  -10.  We  see  from  Figs.  10  and  12  that,  although  the  value  Tq  is  unstable,  the 
estimate  T(t)  tends  to  undergo  a  self-excited  oscillation  about  the  equilibrium  value 
indicated  by  the  dotted  lines.  The  evidence  favors  the  conjecture  that  in  these  cases  the 
nonlinear  equation  (2.11)  may  exhibit  a  Hopf  type  bifurcation  as  the  parameter  Tq  passes 
from  .1  (10  Hz.)  to  .0667  (15  Hz.).  Detailed  analysis  of  this  possibility  must  await 
later  treatment. 

Figures  13  through  16  show  experience  in  the  15  Hz  case  with  k^  and  k2  selected 
so  that  X  =  -14  (Figs.  13  and  14)  and  so  that  X  =  -20  (Figs.  15  and  16).  We  see  that 


the  performance  improves  as  (5.4)  is  made  progressively  more  stable,  in  agreement  with  the 
conjectures  of  Section  3.  The  small  residual  oscillation  evident  in  Fig.  15  is  probably 
due  to  the  fact  that  PERIOD  does  not  provide  an  exact  estimate  even  when  the  corresponding 
continuous  process  associated  with  minimization  of  (2.4)  or  (2.6)  is  asymptotically 
stable. 


Figures  17  and  IS  show  variational  AT  solutions,  obtained  in  the  manner  described 


at  the  end  of  Section  3,  for  the  15  Hz  case  with  X  =  -9  and  X  =  -14,  respectively. 
Because  T(t)  does  not  converge  to  Tq  in  the  first  case  (cf.  Fig.  12),  even  the 
variational  solutions  are  not  sinusoidal. 

In  Fig.  18,  corresponding  to  X  -  -14,  T(t)  converges  to  Tq  «  .0667  (cf.  Fig.  14) 
and  the  corresponding  variational  solution  tends  to  zero  in  a  convincing  exponentially 
damped  sinusoidal  manner,  this  behavior  becoming  more  convincing  as  t  increases.  It  is 
of  interest  to  estimate  the  frequencies  and  damping  factors  here  and  compare  them  with  the 
eigenvalues  of  N(aQ),  F(aQ).  Analyzing  the  data  plotted  in  Fig.  18  one  obtains  the 
estimate  T  »  .66  and,  comparing  the  successive  amplitudes,  we  see  that  the  oscillation 


and  its  eigenvalues  may  be  computed  to  be 
-2.77  ±  1105.32,  -10.82  *  15.9  . 

None  of  these  correspond  to  (5.5)  and  we  conclude  that  the  dynamics  exhibited  in  the 
identification  process  arise  from  a  different  source  -  which,  on  the  basis  of  our  earlier 
investigations  in  this  paper,  we  believe  to  be  the  functional  equation  (2.22) 
(equivalently  (2.19),  (2.20)). 

In  Figures  19-26  we  indicate  the  effect  of  varying  the  parameters  y  and  e 


(cf.  (1.9),  (2.4))  while  leaving  the  feedback  parameters  in  K  fixed  at  the  values  which 


j — ■r;  y.1  y 1  •* "  ■* •  r  .'i 


I 

.  produced  Figures  13  and  14  with  Y  =  .9  and  e  “  1000.  This  corresponds  to  the  double 

[  eigenvalue  X  «  -14  for  (5.4).)  Here  it  needs  to  be  explained  that  the  value  "GAMMA" 

referred  to  on  the  figure  heading  corresponds  to  e  where  h  is  the  length  of  the 

sampling  interval  used  by  PERIOD.  For  all  cases  studied  here  h  is  ten  times  the  H 
value  shown  in  the  figure  heading;  thus  h  =  .01  and 

i 

—  .  0  1 Y 

GAMMA  *  .8  3  e  1 
log  .8 

♦  Y  ”  - - —  =  -22.3  , 

-.01 

GAMMA  =  .9  +  Y  =  -10.5 
GAMMA  =  .95  ♦  Y  =  -5.13  . 

As  we  see  by  comparison  of  Figs.  19  and  20  with  13  and  14,  performance  is  substantially 
degraded  by  discounting  past  values  too  much,  corresponding  to  y  =  “5.13  on  the  other 
hand,  gives  substantially  better  performance  than  Y  “  -10.5. 

For  Figs.  23  and  24  we  have  set  GAMMA  =.9  again  but  have  increased  t  to  2000 
rather  than  the  earlier  1000.  The  improvement  over  the  results  in  Figs.  13  and  24  is 
again  quite  marked.  To  obtain  the  results  of  Figs.  25  and  26  we  "pulled  out  all  the 
stops"  and  set  GAMMA  -  .95,  e  =  2000  and  chose  feedback  parameters  corresponding  to  a 
double  eigenvalue  X  -  -20  for  (5.4).  Here  we  finally  obtain  the  sort  of  sinusoidal 
disturbance  rejection  one  would  hope  for. 

It  is  clear  from  these  numerical  studies  that  the  procedure  we  have  described  here 
can  be  effective  for  vibration  suppression  under  certain  circumstances.  We  have  also  made 
a  case  in  this  article  for  the  proposition  that  the  near  equilibrium  plant/coiqpensator/ 
identifier  dynamics  are  governed  by  a  functional  equation  involving  an  infinite  delay  and 
periodic  coefficients.  It  remains  for  other  investigations  to  develop  the  detailed 
relationships  between  the  design  parameters  and  the  dynamics  of  solutions  of  this 
functional  equation  upon  which  a  design  methodology  for  effective  performance  might  be 


based. 


0010  GAMMA  =  0000  FIJI 


3  5KD2=  -0  7  F'RFIQ  = 


tl 


6.  sketch  of  the  pkoof  of  THEOKEM  5. 

Let  z(t),  w( t , s ) ,  etc.,  be  as  in  the  statement  of  Theorem  5.  For  any  t  >  0  we 

have 

t 

z(t)  =  J  W(t,s)z(s)ds  = 


t  kT  0 

J  W(t,s)z(s)ds  +  /  W(t,s)z(s)ds  +  J  W(t,s)z(s)ds  , 
kT  0  -® 

2 

where  k  is  the  largest  integer  such  that  kT  <  t.  We  define  z^  €  L^IO.T)  by 

z£(t)  -  z( ( 4- 1  )T  +  t),  t-C  [0,T],  4  >  1  , 

~  2 

and  we  define  z  0  e  L  [0,t]  by 
“  x  m 

z  4<t)  -  z(-<4+1)T  +  Tj,  t  e  [0,Tl,  4  >  0  . 

Then,  with 

t  =  kT  +  T,  s  =  4T  +  O,  s  e  [4T,(4+1)T]  , 

(6.1)  yields 

T  k  T 

-  /  w(kT  +  T,kT  +  o  )z.  .  (o)do  -  l  j  W(kT  +  t  ,  ( j  -  1  )T  +  o)z. 

k+1  0  K  '  j-1  0  3 


(6.1) 


(o  )do 


=  l  j  W(kT  +  T,  -(4  +  1 )T  +  0)z  . ( 0 ) do  . 
4-0  0 


Using  the  periodicity  relation  (3.7)  we  have 

W(kT  +  T,  (j  -  1 )T  +  O)  =  w(t,  -  (k  -  j )T  +  o)  5  Wk.j(T,0) 
for  any  integer  j  <  k.  Then  (6.2)  becomes 


(6.2) 


k+1 


T  k  T 

(T)  -  /  w0(T,o)zk+1(o)do  -IS  Wk+1_j(T,o)z . (o)do 

0  j=1  0  J  J 


1  J  W^^Oz.^ddo  * 
4=0  0 


(6.3) 
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4 


The  conditions  (3.7),  (3.8)  show  that  the  last  sum  converges  in  L^tO.T).  we  may  write 

2 

(6.3)  as  a  vector  linear  recursion  equation  in  L^tO.T] : 

k  ®  __ 

(1  "  P0,zk+1  '  \  Pk+1-j  zj  “  ^  Pk+1+i  z  l  t6*4) 

j=1  J  J  i«0 


(Pqz)(t)  -  /  W-  (t  ,o  )z(c  )do , 
0 


(P.  z)(t)  -  /  W.  (T,0)z(0)da,  k  »  1,2,3  .... 


As  is  well  known  ([3])  I  -  PQ  is  boundedly  invertible,  Pq,  P1 ,  P2,  ...  are  all 
compact .  Then 

Q_k  -  -(I  -  P0>"\.  *  “  1.2,3,  ...  (6.5) 

are  all  compact  and,  keeping  (3.7)  in  mind,  it  may  be  seen  that  for  some  positive 
number  D 

«Q_kl  <  De“kcT,  k  -  1,2,3 .  (6.6) 

Then  (6.4)  can  be  written, with  an  obvious  re-indexing,  as 
-1  -k 

Zk  +  1  QiZt+k  +  ^  24Zl+k  '  k  ~  1»2»3»  ....  (6.7) 

k  l— k+1  k  i— ®  1  *  k 

Given  a  sequence  {y ^ | —  ®  <  k  <  *}  C  H,  where  H  is  a  Hilbert  space,  and  supposing 

that 

lykl  <  M+(T+)k,  k  -  1,2,3,  ...»  (6.8) 

^  M  (Tf“)  .  k  “  0,-1  ,-2,-3, . . . ,  (6.9) 

where  M+,M  ,y+,Y  are  all  positive  numbers  and  Y+  >  Y  ,  we  define  the  bilateral 


"Z-transform"  (discrete  Laplace  transform)  of  {yk}  by 
00 

n<Y>  -  l  y.X-k  =  h+(X),  |Xj  >  y+  , 


).  y  kX  =~T1  (X),  I X  j  <  Y~ 


(6.10) 


Clearly,  n(X)  is  analytic  in  neighborhoods  of  both  0  and  00 .  In  many  cases 
n  (1)  and  n-( X )  are  analytic  continuations  of  each  other.  For  example,  if  for  all 


integers  k 

yk  =  uky0»  y0e  h'  v  *  0  e 

then  with  Y+  =  Y  =  | V \ ,  M+  =  M  =1,  all  of  the  above  are  valid  for 


n+(X)  =  n(X),  111  >  In|,  n— ( X )  =  n(X),  j X 1  <  |u|,  n(X>  =  1 


X  -  u  * 


If,  correspondingly,  i  ~  <  k  <  °®|  ,  is  a  sequence  of  bounded  operators  on  H 


such  that 


IQk«  <  B+(P+)"k,  k  =  1,2,3 


lQ_kl  <  B  (p  )  ,  k  =  0,-1, -2,  ..., 


(6.11) 

(6.12) 


where  B+,  B-,  p+,  p  are  all  positive  (and  p+  >  p  -  <  Y  in  our  application),  we 


may 


define  the  discrete  Fourier  transform  of  { Q.  }  by 


oa>  =  1  • 

k“-“ 


(6.13) 


Clearly  the  series  converges  and  Q(X)  is  a  holomorphic  operator  valued  function  for 
p  <  |X|  <  p+  .  if  0^  =  0  for  all  positive  k,  then  P+  may  be  taken  to  be  ■  and 
Q(X)  will  be  holomorphic  for  |X|  >  p  ,  including  X  = 

The  convolution  of  {Qk}  =  Q  and  {y^}  =  y  is  defined  by 


5  u'y’i  =  J  Vk+f* 

k=-“ 


(6.14) 


the  sum  being  convergent  when  (6.8),  (6.9),  (6.11)  and  (6.12)  apply  and 
p+  >  Y+,  P  -  <  Y  ,  as  we  suppose.  To  anyone  familiar  with  transforms  of  this  type  the 
first  question  occurring  concerns  the  relationship  of  the  transform  $(X)  of  {f^}  to 
the  transforms  Q(X),  n(X).  The  answer  is  easy  but  not  completely  obvious.  Let 
r+  and  T  be  positively  oriented  circles  centered  at  X  =  0  with  radii  r+,  r“, 

Y+  <  r+  <  p+ ,  p  <  r  <  y  •  Then,  as  we  show  in  the  more  complete  discussion  [7),  if 


r  =  r  -  r  and  X  lies  in  the  exterior  of  the  annular  regions  between  T  and  T  , 


that  meets  no  singularities  of  Q(X).  Then,  applying  (6.16)  to  the  zj,  of  (6.7)  we 


zh  =  2ii  i  n(X,xk-1dX  +31-  /  n(X,xk-1cx 
6  1 6 


c,F  +  zk,r"+6  '  k  “1'2'3 . 


where  =  T  -  .  From  (6.19)  and  (6.21)  it  is  clear  that 

lz^,r  +  Si  <  M(r  +  6)k 

where  M  is  a  constant  which  may  be  bounded  in  terms  of  h,  hence  in  terms  of  the 
z  ,  A  =0,1,2,  ....  On  the  other  hand 

y  v 

z  =.  t:  rX.  Res  r\  ( X  . )  . 
k,F  X2e  Intrfi  3  j 

In  the  case  of  a  simple  eigenvalue  X..  with  one  dimensional  null  space,  which  is  all  we 
will  study  here,  from  the  formula  (6.19)  for  n  it  may  be  seen  that 


(♦.,*Q(C)MC)d* 


Res  MX.)  -  (4,.,Q*(x.)^.)  J- 


(6.22) 


where  is  a  non-zero  vector  in  the  one  dimensional  null  space  of  Q(Xj)  and  <|k  is  a 

corresponding  vector  in  the  null  space  of  Q(X.)*  such  that  (\Jj  .  ,$  .  )  _  =  0.  We  see 

3  33  l/[0,T] 

in  any  case  that  z  is  a  sum  of  the  form 
kfF 


z,_  „  -  ,  _  T — .  r  Xk  c.<p.  ,  k  =  1,2,3,  .... 
k , F  x_.e  Int  J  ]'j 


(6.23) 


The  corresponding  solutions  zF(t),  ( t )  of  (3.6)  (or  (6.1))  are  obtained  by  inverting 

-6t 

the  transformations  which  follow  (6.1).  The  term  e  of  (3.11)  is  identified  with 

r  +6.  It  is  greater  than  e_c^  which  is  identified  with  r”.  Thus  z„(t)  satisfies 

P 

(3.11). 
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Since  *p(t)  is  a  solution  of  (6.2),  the  form  of  that  equation  shows  that  zp(t) 
must  be  a  continuous  function.  The  form  (6.23)  then  implies  that  zp  on  any  interval 
[kT,(k  +  1 )T]  is  Xj  times  the  corresponding  value  of  zp  on  [(k  -  1)T,kT].  From  this 
is  clear  that 

(t)  -  Xj  *  ,0)  . 

At 

We  identify  e  in  (3.12)  with  X^  and  P(T)  with  the  T  -  periodic  extension  of 

e  ^^(t)  (which  satisfies 

—At  ->1  — x  o 

P(T)  -  e  ^(T)  -  X.  (T)  -  ifr^O)  -  e  *.(0>  -  P(0)) 

Thus,  modulo  the  usual  remarks  which  most  apply  to  non-simple  poles  of  Q(X)  \  which 
lead  to  solutions  of  the  form  (3.13),  we  have  completed  the  proof  of  Theorem  5.  Further 
details  may  be  found  in  [7] .  The  main  point  of  the  theorem  is  that  the  dominant  solutions 
of  (3.6)  (or  (6.1))  are  those  associated  with  the  larger  singularities  of  Q(X)  and  those 
solutions  are  of  the  type  (3.12)  or  (3.13). 
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